The Molecular Transducers of Physical Activity Consortium (MoTrPAC) conducted the first comprehensive, organism-wide, multi-omic investigation of endurance exercise training effects in rats. This large-scale study involved six-month-old male and female Fischer 344 rats subjected to progressive treadmill endurance exercise training over 1, 2, 4, or 8 weeks. Tissues were collected from these rats 48 hours after the last exercise bout to capture the adaptations resulting from the training, with sex-matched sedentary, untrained rats serving as controls. Utilizing multiple ‘omic’ platforms including transcriptomics, epigenomics, proteomics, and metabolomics—across a wide range of tissues, the study provides a rich dataset that captures the temporal effects of endurance exercise training at multiple biological levels. This dataset offers invaluable insights into the molecular mechanisms underlying exercise-induced adaptations. Find the open access manuscript published on Nature.
To facilitate access to and analysis of this extensive dataset, two essential R packages have been developed:
In this tutorial, we will explore how to utilize the
MotrpacRatTraining6moData package to access and navigate
the extensive dataset generated by the endurance exercise study, and
employ the MotrpacRatTraining6mo package to perform
analyses and gain insights into the molecular responses to endurance
training in young rats.
Note: This notebook is based on the vignettes available in MotrpacRatTraining6moData and MotrpacRatTraining6mo vignette.
Install both packages with devtools:
if (!require("devtools", quietly = TRUE)){
install.packages("devtools")
}
options(timeout=1e5) # extend the timeout for downloading large files
# First the data package, which might take a long time
devtools::install_github("MoTrPAC/MotrpacRatTraining6moData")
# Then the analysis package:
devtools::install_github("MoTrPAC/MotrpacRatTraining6mo")Once we install this package, we can load these and other libraries used in this tutorial:
library(MotrpacRatTraining6moData)
library(MotrpacRatTraining6mo)
library(ggplot2) # for plots in this tutorial
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, unionlibrary(impute) # For KNN imputation
library(tidyverse)
#> ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
#> ✔ forcats 1.0.0 ✔ stringr 1.5.1
#> ✔ lubridate 1.9.3 ✔ tibble 3.2.1
#> ✔ purrr 1.0.2 ✔ tidyr 1.3.1
#> ✔ readr 2.1.5
#> ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
#> ✖ dplyr::filter() masks stats::filter()
#> ✖ dplyr::lag() masks stats::lag()
#> ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsDetails of the experimental design can be found in the supplementary methods of the Nature paper.
Briefly, 6-month-old young adult rats were subjected to progressive endurance exercise training for 1, 2, 4, or 8 weeks, with tissues collected 48 hours after the last training bout. Sex-matched sedentary, untrained rats were used as controls. Whole blood, plasma, and 18 solid tissues were analyzed using genomics, proteomics, metabolomics, and protein immunoassay technologies, with most assays performed in a subset of these tissues. Depending on the assay, between 3 and 6 replicates per sex per time point were analyzed.
The PHENO object represents phenotypic data from the
MoTrPAC endurance exercise training study in six-month-old rats. It is a
comprehensive dataset containing 5,955 rows and 510
variables, with each row corresponding to a unique sample
identified by viallabel. The dataset captures detailed
information on the animals, their training regimens, specimen
collection, and various physiological metrics. Here is a summary of the
major components:
Identifiers: Unique identifiers such as
pid, bid, viallabel, and
labelid are used to link samples to individual animals and
their specimen labels.
Animal Information: Variables such as
registration.d_birth, registration.sex, and
registration.weight provide basic details about the rats,
including birth date, sex, and weight.
Registration and Housing: Data about the
arrival, housing conditions, registration
(registration.d_arrive,
registration.cagenumber), and light conditions
(registration.d_reverselight) are included to document the
environmental conditions in which animals were kept.
Intervention and Randomization: Information
about intervention groups (e.g., control or training) is captured in
variables like key.intervention,
key.anirandgroup, and key.protocol, indicating
how each rat was assigned to different treatment conditions.
Familiarization and Training Data: Detailed
records of treadmill familiarization
(familiarization.d_treadmillbegin) and progressive
endurance training are included. Training data cover up to 40
days, with variables for each day’s treadmill speed, incline,
and exercise time (training.dayX_treadmillspeed,
training.dayX_timeontreadmill).
VO2 Max and NMR Testing: Information from
VO2 max tests and NMR body composition
tests is included. Variables such as
vo2.max.test.vo2_max and nmr.testing.nmr_fat
capture maximum oxygen uptake, percent body fat, and other key metabolic
measures.
Specimen Collection: Specimen collection details
are provided, including dates
(specimen.collection.d_visit), times of collection
(specimen.processing.t_collection), and types of tissues
collected
(specimen.processing.sampletypedescription).
Terminal Measures: Terminal metrics, such as
body weight at sacrifice and tissue weights
(terminal.weight.bw, terminal.weight.sol),
document the final physical characteristics of the animals after
completion of the study.
Calculated Variables: Derived variables, such as
changes in body composition
(calculated.variables.pct_body_fat_change), lactate levels,
and vo2_max, provide insights into physiological
adaptations resulting from the exercise intervention.
Custom Variables: Variables like
sacrificetime, intervention, and
time_to_freeze have been added for simplified grouping and
analysis, reflecting key time points, interventions, and
time-to-processing of samples.
pid: Participant Identifier. A unique
8-digit identifier assigned to each animal subject (rats in this study).
All samples (viallabel) coming from the same animal have
the same pid. This variable is crucial when combining
results from multiple assays with phenotypic data, allowing for
animal-specific longitudinal analysis.viallabel: Vial Label ID. A unique
11-digit code assigned to each sample vial. This ID is present across
all related results and metadata files, serving as a key to link the
quantitative results with the phenotypic data.sex: Sex of the animal. Represented as
"female" for Female and "male". This variable
is critical for any analysis that aims to determine sex-based
differences.group: A simplify version of
study_group_timepoint that only includes the group and time
points: control (sedentary), and 1w,
2w, 4w, and 8w (exercise)tissue: tissue name (short)In summary, the PHENO object provides a comprehensive overview of each animal’s involvement in the study, including the conditions under which they were kept, the specific training they underwent, and the physiological changes observed as a result of the endurance training intervention. This dataset allows for in-depth analysis of the effects of physical activity at various biological levels.
First, let’s select the ids of the rats used in this study:
# Load the proteomics data, for example:
data_prot <- MotrpacRatTraining6mo::combine_normalized_data(
assays = c("PROT","UBIQ","PHOSPHO","ACETYL"),
exclude_outliers = TRUE, )
#> PROT_CORTEX_NORM_DATA
#> PROT_HEART_NORM_DATA
#> PROT_KIDNEY_NORM_DATA
#> PROT_LIVER_NORM_DATA
#> PROT_LUNG_NORM_DATA
#> PROT_SKMGN_NORM_DATA
#> PROT_WATSC_NORM_DATA
#> UBIQ_HEART_NORM_DATA
#> UBIQ_LIVER_NORM_DATA
#> PHOSPHO_CORTEX_NORM_DATA
#> PHOSPHO_HEART_NORM_DATA
#> PHOSPHO_KIDNEY_NORM_DATA
#> PHOSPHO_LIVER_NORM_DATA
#> PHOSPHO_LUNG_NORM_DATA
#> PHOSPHO_SKMGN_NORM_DATA
#> PHOSPHO_WATSC_NORM_DATA
#> ACETYL_HEART_NORM_DATA
#> ACETYL_LIVER_NORM_DATA
# Get the unique pid (rat id) values from the proteomics data
rat_omics <- setdiff(names(data_prot), c("feature", "feature_ID", "tissue", "assay"))
# Count unique 'pid' by 'group' and 'sex'
count_data <- PHENO %>%
group_by(group, sex) %>%
summarise(unique_pid_count = n_distinct(pid),
.groups = 'drop') %>%
ungroup()
# Create the grouped bar chart
p <- ggplot(count_data, aes(x = group, y = unique_pid_count, fill = sex)) +
geom_bar(stat = "identity", position = position_dodge()) +
labs(title = "Animal count",
subtitle = "All animals processed at the Animal Site",
x = "Study Group Time Point",
y = "Count of Unique Animals (pid)",
fill = "Sex") +
theme_linedraw() +
coord_flip()
print(p)
# Animals used in this study
# Filter by pid (i.e., the animal id)
count_data <- PHENO %>%
filter(pid %in% rat_omics) %>%
group_by(group, sex) %>%
summarise(unique_pid_count = n_distinct(pid),
.groups = 'drop') %>%
ungroup()
# Create the grouped bar chart
p2 <- ggplot(count_data, aes(x = group, y = unique_pid_count, fill = sex)) +
geom_bar(stat = "identity", position = position_dodge()) +
labs(title = "Animal Count",
subtitle = "Only animals used for molecular profile",
x = "Study Group Time Point",
y = "Count of Unique Animals (pid)",
fill = "Sex") +
theme_linedraw() +
coord_flip()
print(p2)# Count total number of unique 'viallabel' by 'tissue'
viallabel_data <- PHENO %>%
filter(pid %in% rat_omics) %>% # only include animals used in the omics data
group_by(tissue, group) %>%
summarise(total_viallabel_count = n_distinct(viallabel), .groups = 'drop') %>%
ungroup()
p <- ggplot(viallabel_data, aes(x = tissue, y = total_viallabel_count, fill = group)) +
geom_bar(stat = "identity") +
labs(title = "Total Number of Samples collected for Each Tissue",
x = "Tissue",
y = "Total Count of Vial Labels",
fill = "Group") +
theme_minimal() +
theme(legend.position = "right") + coord_flip()
print(p)# Select variables from PHENO that start with "training" and end with "_weight"
weight_cols <- grep("^training.*_weight$", names(PHENO), value = TRUE)
# Select those columns along with "pid", "sex", and "group"
selected_cols <- c("pid", "sex", "group", weight_cols)
# Get unique values
unique_pheno <- unique(PHENO[, selected_cols])
# Reshape data to long format for plotting
pheno_long <- unique_pheno %>%
dplyr::filter(pid %in% rat_omics) %>%
tidyr::pivot_longer(cols = starts_with("training"), names_to = "training_session", values_to = "weight")
# Define the desired order of training sessions
training_session_order <- c("training.day1_weight",
"training.day6_weight",
"training.day11_weight",
"training.day16_weight",
"training.day21_weight",
"training.day26_weight",
"training.day31_weight",
"training.day36_weight")
# Convert training_session to a factor with the specified order
pheno_long$training_session <- factor(pheno_long$training_session,
levels = training_session_order)
# Plot the jittered plot for each training session, faceted by sex and colored by group
p1 <- ggplot(pheno_long, aes(x = training_session, y = weight, color = group)) +
geom_jitter(width = 0.2, size = 2, alpha = 0.6, na.rm = TRUE) +
facet_wrap(~sex) +
theme_bw() +
labs(title = "Weight by Training Session",
x = "Training Session",
y = "Weight") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
print(p1)
p2 <- ggplot(pheno_long, aes(x = training_session, y = weight, color = group)) +
geom_boxplot(alpha = 0.3, outlier.shape = NA, na.rm = TRUE, linewidth = 0.5) +
facet_wrap(~sex, scales = "free") +
theme_bw() +
labs(title = "Weight by Training Session",
x = "Training Session",
y = "Weight") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
print(p2)
p3 <- ggplot(pheno_long, aes(x = training_session, y = weight, group = group, color = group, fill = group)) +
stat_summary(fun = median, geom = "point", size = 4, shape = 21, alpha = 0.7, na.rm = TRUE) +
stat_summary(fun = median, geom = "line", size = 1, alpha = 0.7, na.rm = TRUE) +
facet_wrap(~sex, scales = "free") +
theme_linedraw() +
labs(title = "Mean Weight by Training Session",
x = "Training Session",
y = "Mean Weight") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.VO2max measurements were taken at the start for all time points, and at the end for the 4-week, 8-week, and control training sessions.
# colnames(PHENO)[grep("vo2*d_visit", colnames(PHENO))]
# Select those columns along with "pid", "sex", and "group"
selected_cols <- c("pid",
"sex",
"group",
"vo2.max.test.vo2_max_1",
"vo2.max.test.vo2_max_2")
# Get unique values
unique_pheno <- unique(PHENO[, selected_cols])
# Reshape data to long format for plotting
pheno_long <- unique_pheno %>%
dplyr::filter(pid %in% rat_omics) %>%
tidyr::pivot_longer(cols = starts_with("vo2"),
names_to = "training_session",
values_to = "vo2max")
p1 <- ggplot(pheno_long,
aes(x = training_session, y = vo2max, color = group)) +
geom_boxplot(alpha = 0.3, outlier.shape = NA, na.rm = TRUE) +
# geom_jitter(width = 0.2, size = 2, alpha = 0.6, na.rm = TRUE) +
facet_wrap(~sex, scales = "free") +
theme_bw() +
labs(title = "VO2max distributions at first and last training session",
x = "Training Session",
y = "VO2max") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
print(p1)
p2 <- ggplot(pheno_long, aes(x = training_session, y = vo2max, group = group, color = group, fill = group)) +
stat_summary(fun = mean, geom = "point", size = 4, shape = 21, alpha = 0.7, na.rm = TRUE) +
stat_summary(fun = mean, geom = "line", size = 1, alpha = 0.7, na.rm = TRUE) +
facet_wrap(~sex) +
theme_linedraw() +
labs(title = "Mean VO2max at first and last training session",
x = "Training Session",
y = "Mean VO2max") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
print(p2)# Load required packages
library(lubridate)
# Filter the data for training rats
PHENO_training <- PHENO %>%
filter(intervention == "training")
# Pivot the date columns
PHENO_long_date <- PHENO_training %>%
select(pid, group, matches("training\\.day\\d+date$")) %>%
pivot_longer(
cols = -c(pid, group),
names_to = "day",
names_pattern = "training\\.day(\\d+)date$",
values_to = "date"
) %>%
mutate(
day = as.integer(day),
date = dmy(date)
)
# Pivot the treadmillspeed columns
PHENO_long_speed <- PHENO_training %>%
select(pid, group, matches("training\\.day\\d+_treadmillspeed$")) %>%
pivot_longer(
cols = -c(pid, group),
names_to = "day",
names_pattern = "training\\.day(\\d+)_treadmillspeed$",
values_to = "treadmillspeed"
) %>%
mutate(
day = as.integer(day),
treadmillspeed = as.numeric(treadmillspeed)
)
# Merge the two datasets
PHENO_long <- suppressWarnings(left_join(PHENO_long_date, PHENO_long_speed, by = c("pid", "group", "day")) %>% distinct())
# Remove rows where both 'date' and 'treadmillspeed' are NA
PHENO_long <- PHENO_long %>%
filter(!(is.na(date) & is.na(treadmillspeed)))
# Create a new Id combining group and pid
PHENO_long <- PHENO_long %>%
mutate(group_pid = paste0(group, "_", pid))
# # Identify duplicates
# duplicates <- PHENO_long %>%
# select(group_pid, date, treadmillspeed) %>%
# group_by(group_pid, date) %>%
# filter(n() > 1)
#
# # View duplicates
# print(duplicates)
# Reshape the data to wide format
PHENO_wide <- PHENO_long %>%
pivot_wider(
names_from = date,
values_from = treadmillspeed
)
# Aggregate duplicates by taking the mean of treadmillspeed
PHENO_long <- PHENO_long %>%
group_by(group_pid, date) %>%
summarise(
treadmillspeed = mean(treadmillspeed, na.rm = TRUE),
.groups = 'drop'
)
# Reshape the data to wide format
PHENO_wide <- PHENO_long %>%
pivot_wider(
names_from = date,
values_from = treadmillspeed
)
# Convert the data frame to a matrix
PHENO_matrix <- PHENO_wide %>%
column_to_rownames(var = "group_pid") %>%
as.matrix()
# Replace NA values with zeros
PHENO_matrix[is.na(PHENO_matrix)] <- 0
library(pheatmap)
library(viridis)
#> Loading required package: viridisLite
# Create the heatmap
pheatmap(
PHENO_matrix,
cluster_rows = FALSE,
cluster_cols = FALSE,
show_rownames = TRUE,
show_colnames = TRUE,
scale = "none",
fontsize_row = 4,
fontsize_col = 4,
angle_col = 45,
color = viridis::viridis(100),
main = "Treadmill Speed Over Time for Training Rats"
)Tip: To learn more about any data object, use
? to retrieve the documentation, e.g.,
?METAB_FEATURE_ID_MAP. Note that
MotrpacRatTraining6moData must be installed and loaded with
library() for this to work.
It is important to be aware of the tissue and assay abbreviations
because they are used to name many data objects. The vectors of
abbreviations are also available in TISSUE_ABBREV and
ASSAY_ABBREV.
Here is a brief summary of the kinds of data included in the data package:
PHENOFEATURE_TO_GENE, RAT_TO_HUMAN_GENEMETAB_FEATURE_ID_MAP, METHYL_FEATURE_ANNOT
(GCP only), ATAC_FEATURE_ANNOT (GCP only),
PROT_FEATURE_ANNOT (GCP only),
PHOSPHO_FEATURE_ANNOT (GCP only),
UBIQ_FEATURE_ANNOT (GCP only),
ACETYL_FEATURE_ANNOT (GCP only),
TRNSCRPT_FEATURE_ANNOT (GCP only)TRNSCRPT_META, ATAC_META,
METHYL_META, IMMUNO_META,
PHOSPHO_META, PROT_META,
ACETYL_META, UBIQ_METATRNSCRPT_LIVER_RAW_COUNTS. Note that
epigenetic data (ATAC and METHYL) must be downloaded from Google Cloud
Storage. See more details here.TRNSCRPT_SKMGN_NORM_DATAHEART_PROT_DAOUTLIERSTRAINING_REGULATED_FEATURESGRAPH_PW_ENRICHA list of all the available data objects and a brief description are available here.
# List all data objects in the package
all_objects <- data(package = "MotrpacRatTraining6moData")$results[, "Item"]
# Filter the objects that contain the string "NORM" and exclude specific unwanted objects
exclude_objects <- c("IMMUNO_NORM_DATA_FLAT", "IMMUNO_NORM_DATA_NESTED",
"METAB_NORM_DATA_FLAT", "METAB_NORM_DATA_NESTED",
"TRAINING_REGULATED_NORM_DATA", "TRAINING_REGULATED_NORM_DATA_NO_OUTLIERS")
norm_objects <- all_objects[grepl("NORM_DATA", all_objects) & !(all_objects %in% exclude_objects)]Use load_sample_data() to load sample-level data for a
specific ome and tissue. Here we fetch various forms of the sample-level
RNA-Seq (TRNSCRPT) data for brown adipose tissue (BAT) as an
example.
# Load RNA-seq raw counts for brown adipose tissue
data <- load_sample_data("BAT", "TRNSCRPT", normalized = FALSE)
# Load the normalized RNA-seq data for brown adipose tissue instead
data <- load_sample_data("BAT", "TRNSCRPT")
# Load the normalized RNA-seq data for brown adipose tissue, but exclude sample outliers
data <- load_sample_data("BAT", "TRNSCRPT", exclude_outliers = TRUE)
# Load the normalized RNA-seq data for brown adipose tissue for training-regulated features only
data <- load_sample_data("BAT", "TRNSCRPT", training_regulated_only = TRUE)load_sample_data() will tell you if the specified
dataset doesn’t exist.
data = load_sample_data("VENACV", "PROT")
#> Warning in load_sample_data("VENACV", "PROT"): No data returned for tissue
#> VENACV and assay PROT with current arguments.load_sample_data() will also download epigenetic data
from Google Cloud Storage.
Note: Epigenetic data require substantially more memory than other omes.
# Load ATAC-seq raw counts for hippocampus, excluding outliers
data = load_sample_data("HIPPOC",
"ATAC",
exclude_outliers = TRUE,
normalized = FALSE,
scratchdir = "/tmp")We can plot the normalized sample-level data for a single feature
using plot_feature_normalized_data(). All of the following
examples are different ways to plot the same feature.
# Two different ways to plot the same results for one feature
plot_feature_normalized_data(feature = "PROT;SKM-GN;NP_786937.1",
exclude_outliers = TRUE,
scale_x_by_time = FALSE)
plot_feature_normalized_data(assay = "PROT",
tissue = "SKM-GN",
feature_ID = "NP_786937.1",
exclude_outliers = TRUE,
scale_x_by_time = FALSE)
plot_feature_normalized_data(assay = "PROT",
tissue = "SKM-GN",
feature_ID = "NP_786937.1",
facet_by_sex = TRUE)combine_normalized_data() is a wrapper for
load_sample_data() that returns combined sample-level
normalized data for multiple datasets. Note that the sample-specific
vial labels used as column names for most of the sample-level data are
replaced with rat-specific participant IDs (PIDs) to allow measurements
from multiple datasets for the same animal to be concatenated.
# Return all normalized RNA-seq data
data_tra <- combine_normalized_data(assays = "TRNSCRPT")
# Return all normalized proteomics data. Exclude outliers
data_prot <- combine_normalized_data(assays = c("PROT","UBIQ","PHOSPHO","ACETYL"),
exclude_outliers = TRUE)
ratid_tra <- setdiff(names(data_tra), c("feature", "feature_ID", "tissue", "assay"))
ratid_prot <- setdiff(names(data_prot), c("feature", "feature_ID", "tissue", "assay"))
# Return normalized ATAC-seq data for training-regulated features
data <- combine_normalized_data(assays = "ATAC",
training_regulated_only = TRUE)
# Return all non-epigenetic data
# Note that the "include_epigen" argument is FALSE by default
data <- combine_normalized_data()Similarly, combine_da_results() concatenates
differential analysis results for multiple datasets.
# Return all non-epigenetic differential analysis results,
# including meta-regression results for metabolomics
res <- combine_da_results()
# Return METHYL and ATAC differential analysis results for gastrocnemius
# res = combine_da_results(tissues="SKM-GN",
# assays=c("ATAC","METHYL"),
# include_epigen=TRUE)Finally, list_available_data() returns a list of all of
the available data objects in the specified package.
list_available_data("MotrpacRatTraining6moData")
#> [1] "ACETYL_HEART_DA"
#> [2] "ACETYL_HEART_NORM_DATA"
#> [3] "ACETYL_LIVER_DA"
#> [4] "ACETYL_LIVER_NORM_DATA"
#> [5] "ACETYL_META"
#> [6] "ASSAY_ABBREV"
#> [7] "ASSAY_ABBREV_TO_CODE"
#> [8] "ASSAY_CODE_TO_ABBREV"
#> [9] "ASSAY_COLORS"
#> [10] "ASSAY_ORDER"
#> [11] "ATAC_BAT_NORM_DATA_05FDR"
#> [12] "ATAC_HEART_NORM_DATA_05FDR"
#> [13] "ATAC_HIPPOC_NORM_DATA_05FDR"
#> [14] "ATAC_KIDNEY_NORM_DATA_05FDR"
#> [15] "ATAC_LIVER_NORM_DATA_05FDR"
#> [16] "ATAC_LUNG_NORM_DATA_05FDR"
#> [17] "ATAC_META"
#> [18] "ATAC_SKMGN_NORM_DATA_05FDR"
#> [19] "ATAC_WATSC_NORM_DATA_05FDR"
#> [20] "FEATURE_TO_GENE"
#> [21] "FEATURE_TO_GENE_FILT"
#> [22] "GENE_UNIVERSES"
#> [23] "GRAPH_COMPONENTS"
#> [24] "GRAPH_PW_ENRICH"
#> [25] "GRAPH_STATES"
#> [26] "GROUP_COLORS"
#> [27] "IMMUNO_ADRNL_DA"
#> [28] "IMMUNO_BAT_DA"
#> [29] "IMMUNO_COLON_DA"
#> [30] "IMMUNO_CORTEX_DA"
#> [31] "IMMUNO_HEART_DA"
#> [32] "IMMUNO_HIPPOC_DA"
#> [33] "IMMUNO_KIDNEY_DA"
#> [34] "IMMUNO_LIVER_DA"
#> [35] "IMMUNO_LUNG_DA"
#> [36] "IMMUNO_META"
#> [37] "IMMUNO_NORM_DATA_FLAT"
#> [38] "IMMUNO_NORM_DATA_NESTED"
#> [39] "IMMUNO_OVARY_DA"
#> [40] "IMMUNO_PLASMA_DA"
#> [41] "IMMUNO_SKMGN_DA"
#> [42] "IMMUNO_SKMVL_DA"
#> [43] "IMMUNO_SMLINT_DA"
#> [44] "IMMUNO_SPLEEN_DA"
#> [45] "IMMUNO_TESTES_DA"
#> [46] "IMMUNO_WATSC_DA"
#> [47] "METAB_ADRNL_DA"
#> [48] "METAB_ADRNL_DA_METAREG"
#> [49] "METAB_BAT_DA"
#> [50] "METAB_BAT_DA_METAREG"
#> [51] "METAB_COLON_DA"
#> [52] "METAB_COLON_DA_METAREG"
#> [53] "METAB_CORTEX_DA"
#> [54] "METAB_CORTEX_DA_METAREG"
#> [55] "METAB_FEATURE_ID_MAP"
#> [56] "METAB_HEART_DA"
#> [57] "METAB_HEART_DA_METAREG"
#> [58] "METAB_HIPPOC_DA"
#> [59] "METAB_HIPPOC_DA_METAREG"
#> [60] "METAB_HYPOTH_DA"
#> [61] "METAB_HYPOTH_DA_METAREG"
#> [62] "METAB_KIDNEY_DA"
#> [63] "METAB_KIDNEY_DA_METAREG"
#> [64] "METAB_LIVER_DA"
#> [65] "METAB_LIVER_DA_METAREG"
#> [66] "METAB_LUNG_DA"
#> [67] "METAB_LUNG_DA_METAREG"
#> [68] "METAB_NORM_DATA_FLAT"
#> [69] "METAB_NORM_DATA_NESTED"
#> [70] "METAB_OVARY_DA"
#> [71] "METAB_OVARY_DA_METAREG"
#> [72] "METAB_PLASMA_DA"
#> [73] "METAB_PLASMA_DA_METAREG"
#> [74] "METAB_SKMGN_DA"
#> [75] "METAB_SKMGN_DA_METAREG"
#> [76] "METAB_SKMVL_DA"
#> [77] "METAB_SKMVL_DA_METAREG"
#> [78] "METAB_SMLINT_DA"
#> [79] "METAB_SMLINT_DA_METAREG"
#> [80] "METAB_SPLEEN_DA"
#> [81] "METAB_SPLEEN_DA_METAREG"
#> [82] "METAB_TESTES_DA"
#> [83] "METAB_TESTES_DA_METAREG"
#> [84] "METAB_VENACV_DA"
#> [85] "METAB_VENACV_DA_METAREG"
#> [86] "METAB_WATSC_DA"
#> [87] "METAB_WATSC_DA_METAREG"
#> [88] "METHYL_BAT_NORM_DATA_05FDR"
#> [89] "METHYL_HEART_NORM_DATA_05FDR"
#> [90] "METHYL_HIPPOC_NORM_DATA_05FDR"
#> [91] "METHYL_KIDNEY_NORM_DATA_05FDR"
#> [92] "METHYL_LIVER_NORM_DATA_05FDR"
#> [93] "METHYL_LUNG_NORM_DATA_05FDR"
#> [94] "METHYL_META"
#> [95] "METHYL_SKMGN_NORM_DATA_05FDR"
#> [96] "METHYL_WATSC_NORM_DATA_05FDR"
#> [97] "OUTLIERS"
#> [98] "PATHWAY_PARENTS"
#> [99] "PHENO"
#> [100] "PHOSPHO_CORTEX_DA"
#> [101] "PHOSPHO_CORTEX_NORM_DATA"
#> [102] "PHOSPHO_HEART_DA"
#> [103] "PHOSPHO_HEART_NORM_DATA"
#> [104] "PHOSPHO_KIDNEY_DA"
#> [105] "PHOSPHO_KIDNEY_NORM_DATA"
#> [106] "PHOSPHO_LIVER_DA"
#> [107] "PHOSPHO_LIVER_NORM_DATA"
#> [108] "PHOSPHO_LUNG_DA"
#> [109] "PHOSPHO_LUNG_NORM_DATA"
#> [110] "PHOSPHO_META"
#> [111] "PHOSPHO_SKMGN_DA"
#> [112] "PHOSPHO_SKMGN_NORM_DATA"
#> [113] "PHOSPHO_WATSC_DA"
#> [114] "PHOSPHO_WATSC_NORM_DATA"
#> [115] "PROT_CORTEX_DA"
#> [116] "PROT_CORTEX_NORM_DATA"
#> [117] "PROT_HEART_DA"
#> [118] "PROT_HEART_NORM_DATA"
#> [119] "PROT_KIDNEY_DA"
#> [120] "PROT_KIDNEY_NORM_DATA"
#> [121] "PROT_LIVER_DA"
#> [122] "PROT_LIVER_NORM_DATA"
#> [123] "PROT_LUNG_DA"
#> [124] "PROT_LUNG_NORM_DATA"
#> [125] "PROT_META"
#> [126] "PROT_SKMGN_DA"
#> [127] "PROT_SKMGN_NORM_DATA"
#> [128] "PROT_WATSC_DA"
#> [129] "PROT_WATSC_NORM_DATA"
#> [130] "RAT_TO_HUMAN_GENE"
#> [131] "RAT_TO_HUMAN_PHOSPHO"
#> [132] "REPEATED_FEATURES"
#> [133] "REPFDR_INPUTS"
#> [134] "REPFDR_RES"
#> [135] "SEX_COLORS"
#> [136] "TISSUE_ABBREV"
#> [137] "TISSUE_ABBREV_TO_CODE"
#> [138] "TISSUE_CODE_TO_ABBREV"
#> [139] "TISSUE_COLORS"
#> [140] "TISSUE_ORDER"
#> [141] "TRAINING_REGULATED_FEATURES"
#> [142] "TRAINING_REGULATED_NORM_DATA"
#> [143] "TRAINING_REGULATED_NORM_DATA_NO_OUTLIERS"
#> [144] "TRNSCRPT_ADRNL_DA"
#> [145] "TRNSCRPT_ADRNL_NORM_DATA"
#> [146] "TRNSCRPT_ADRNL_RAW_COUNTS"
#> [147] "TRNSCRPT_BAT_DA"
#> [148] "TRNSCRPT_BAT_NORM_DATA"
#> [149] "TRNSCRPT_BAT_RAW_COUNTS"
#> [150] "TRNSCRPT_BLOOD_DA"
#> [151] "TRNSCRPT_BLOOD_NORM_DATA"
#> [152] "TRNSCRPT_BLOOD_RAW_COUNTS"
#> [153] "TRNSCRPT_COLON_DA"
#> [154] "TRNSCRPT_COLON_NORM_DATA"
#> [155] "TRNSCRPT_COLON_RAW_COUNTS"
#> [156] "TRNSCRPT_CORTEX_DA"
#> [157] "TRNSCRPT_CORTEX_NORM_DATA"
#> [158] "TRNSCRPT_CORTEX_RAW_COUNTS"
#> [159] "TRNSCRPT_HEART_DA"
#> [160] "TRNSCRPT_HEART_NORM_DATA"
#> [161] "TRNSCRPT_HEART_RAW_COUNTS"
#> [162] "TRNSCRPT_HIPPOC_DA"
#> [163] "TRNSCRPT_HIPPOC_NORM_DATA"
#> [164] "TRNSCRPT_HIPPOC_RAW_COUNTS"
#> [165] "TRNSCRPT_HYPOTH_DA"
#> [166] "TRNSCRPT_HYPOTH_NORM_DATA"
#> [167] "TRNSCRPT_HYPOTH_RAW_COUNTS"
#> [168] "TRNSCRPT_KIDNEY_DA"
#> [169] "TRNSCRPT_KIDNEY_NORM_DATA"
#> [170] "TRNSCRPT_KIDNEY_RAW_COUNTS"
#> [171] "TRNSCRPT_LIVER_DA"
#> [172] "TRNSCRPT_LIVER_NORM_DATA"
#> [173] "TRNSCRPT_LIVER_RAW_COUNTS"
#> [174] "TRNSCRPT_LUNG_DA"
#> [175] "TRNSCRPT_LUNG_NORM_DATA"
#> [176] "TRNSCRPT_LUNG_RAW_COUNTS"
#> [177] "TRNSCRPT_META"
#> [178] "TRNSCRPT_OVARY_DA"
#> [179] "TRNSCRPT_OVARY_NORM_DATA"
#> [180] "TRNSCRPT_OVARY_RAW_COUNTS"
#> [181] "TRNSCRPT_SKMGN_DA"
#> [182] "TRNSCRPT_SKMGN_NORM_DATA"
#> [183] "TRNSCRPT_SKMGN_RAW_COUNTS"
#> [184] "TRNSCRPT_SKMVL_DA"
#> [185] "TRNSCRPT_SKMVL_NORM_DATA"
#> [186] "TRNSCRPT_SKMVL_RAW_COUNTS"
#> [187] "TRNSCRPT_SMLINT_DA"
#> [188] "TRNSCRPT_SMLINT_NORM_DATA"
#> [189] "TRNSCRPT_SMLINT_RAW_COUNTS"
#> [190] "TRNSCRPT_SPLEEN_DA"
#> [191] "TRNSCRPT_SPLEEN_NORM_DATA"
#> [192] "TRNSCRPT_SPLEEN_RAW_COUNTS"
#> [193] "TRNSCRPT_TESTES_DA"
#> [194] "TRNSCRPT_TESTES_NORM_DATA"
#> [195] "TRNSCRPT_TESTES_RAW_COUNTS"
#> [196] "TRNSCRPT_VENACV_DA"
#> [197] "TRNSCRPT_VENACV_NORM_DATA"
#> [198] "TRNSCRPT_VENACV_RAW_COUNTS"
#> [199] "TRNSCRPT_WATSC_DA"
#> [200] "TRNSCRPT_WATSC_NORM_DATA"
#> [201] "TRNSCRPT_WATSC_RAW_COUNTS"
#> [202] "UBIQ_HEART_DA"
#> [203] "UBIQ_HEART_NORM_DATA"
#> [204] "UBIQ_LIVER_DA"
#> [205] "UBIQ_LIVER_NORM_DATA"
#> [206] "UBIQ_META"If the MotrpacRatTraining6moData
library is attached, you can learn more about any of these data objects
with ?, e.g., ?TISSUE_COLORS,
and you can load data objects into your environment using
data(), e.g.,
data(TRAINING_REGULATED_FEATURES).
Function to perform a PCA on any of the available NORM datasets
# Function to perform PCA and plot the first two PCs
library(impute) # For KNN imputation
do_pca_plot <- function(norm_data_object, object_name = NULL, impute = FALSE, remove_outliers = FALSE) {
# Check if the object name is provided or extract it
if (is.null(object_name)) {
object_name <- deparse(substitute(norm_data_object))
}
# Check if the object name contains 'NORM'
if (!grepl("NORM", object_name)) {
stop("The provided object does not contain 'NORM' in its name. Please provide a valid NORM_DATA object.")
}
# Extract the sample IDs from the list
sample_ids <- setdiff(names(norm_data_object), c("feature", "feature_ID", "tissue", "assay"))
# Subset the data using sample_ids
data_subset <- norm_data_object[, sample_ids]
# Remove columns with more than 80% missing values
valid_cols <- colMeans(is.na(data_subset)) < 0.8
data_subset <- data_subset[, valid_cols]
# Ensure there are enough columns left for PCA
if (ncol(data_subset) < 2) {
stop("Not enough columns remaining after filtering for PCA.")
}
# Impute or clean data based on 'impute' argument
if (impute) {
# Impute missing values using KNN, suppressing messages
sink(file = nullfile())
data_imputed <- impute::impute.knn(as.matrix(data_subset))$data
sink()
} else {
# Remove rows with NA or Inf values
data_imputed <- data_subset[stats::complete.cases(data_subset) & apply(data_subset, 1, function(x) all(is.finite(x))), ]
}
# Ensure there are enough rows left for PCA
if (nrow(data_imputed) < 2) {
stop("Not enough rows remaining after data cleaning for PCA.")
}
# Update the sample_ids based on the final data used in PCA
valid_sample_ids <- colnames(data_imputed)
# Check if sample_ids match 'pid' or 'viallabel' in PHENO and use them accordingly
if (all(valid_sample_ids %in% PHENO$pid)) {
df_sample_ids <- match(valid_sample_ids, PHENO$pid)
} else if (all(valid_sample_ids %in% PHENO$viallabel)) {
df_sample_ids <- match(valid_sample_ids, PHENO$viallabel)
} else {
stop("Sample IDs do not match either 'pid' or 'viallabel' in PHENO.")
}
# Perform PCA on the transposed imputed/cleaned data
pca_result <- stats::prcomp(t(data_imputed))
# Calculate the percentage of variance explained by each PC
variance_explained <- summary(pca_result)$importance[2, ] * 100
# Create a data frame with phenotypic data and the first 3 PCs
df <- data.frame(
group = PHENO[df_sample_ids, "group"],
sex = PHENO[df_sample_ids, "sex"],
pca_result$x[, 1:3] # take the first three principal components
)
# Remove extreme outliers based on IQR for PC1 and PC2 if 'remove_outliers' is TRUE
if (remove_outliers) {
for (pc in c("PC1", "PC2")) {
Q1 <- quantile(df[[pc]], 0.25, na.rm = TRUE)
Q3 <- quantile(df[[pc]], 0.75, na.rm = TRUE)
IQR_value <- Q3 - Q1
# Define the boundaries
lower_bound <- Q1 - 3 * IQR_value
upper_bound <- Q3 + 3 * IQR_value
# Filter the data to remove outliers
df <- df[df[[pc]] >= lower_bound & df[[pc]] <= upper_bound, ]
}
}
# Ensure that there are still data points left to plot
if (nrow(df) < 1) {
stop("No data points left to plot after removing outliers.")
}
# Plot the first two PCs with variance explained in axis labels
p <- ggplot2::ggplot(df, ggplot2::aes(x = PC1, y = PC2, fill = group, shape = sex)) +
ggplot2::geom_point(size = 3, colour = "black") +
ggplot2::scale_fill_manual(values = GROUP_COLORS) +
ggplot2::scale_shape_manual(values = c(male = 21, female = 24)) +
ggplot2::theme_bw() +
ggplot2::guides(fill = ggplot2::guide_legend(override.aes = list(shape = 21))) +
ggplot2::labs(title = paste(object_name),
x = paste0("PC1 (", round(variance_explained[1], 2), "% Variance)"),
y = paste0("PC2 (", round(variance_explained[2], 2), "% Variance)")) +
ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
return(p)
}Perform PCA and plot for all objects in the list of norm objects generated above
for (object_name in norm_objects) {
norm_data_object <- get(object_name)
p <- do_pca_plot(norm_data_object, object_name, impute = FALSE)
print(p)
}The NORM data object for metabolomics contains all
tissues and assays. To perform a PCA, it requires previous filtering and
assay selection. Let’s massage the data to get a PCA for every assay
available for LIVER. Some of the metabolomics assays
require imputation due to a larger number of missing values.
METAB_NORM_ONE_TISSUE <- METAB_NORM_DATA_FLAT %>% filter(tissue == "HEART")
all_assays_one_tissue <- unique(METAB_NORM_ONE_TISSUE$dataset)
for(current_assay in all_assays_one_tissue) {
METAB_NORM_ONE_ASSAY <- METAB_NORM_ONE_TISSUE %>% filter(dataset == current_assay)
METAB_NORM_ONE_ASSAY <- METAB_NORM_ONE_ASSAY %>% select(-dataset)
p <- do_pca_plot(norm_data_object = METAB_NORM_ONE_ASSAY,
object_name = paste0(toupper(current_assay),"_NORM"),
impute = TRUE)
print(p)
}More details about the differential analysis methods are available in the supplementary information of our Nature publication.
Simply put:
# Function to generate a volcano plot for a given data object
create_volcano_plot <- function(data_object) {
# Extract the name of the data object
object_name <- deparse(substitute(data_object))
# Filter out outliers with absolute logFC greater than 10
filtered_data <- data_object %>%
filter(abs(logFC) <= 10)
p <- ggplot(filtered_data, aes(x = logFC, y = -log10(adj_p_value), color = adj_p_value < 0.05)) +
geom_point(alpha = 0.6, size = 1.5) +
scale_color_manual(values = c("TRUE" = "red", "FALSE" = "grey")) +
facet_grid(comparison_group~sex) +
theme_linedraw() +
labs(
title = paste(object_name),
x = "logFC",
y = "-log10(Adjusted P-value)",
color = "Significant (FDR < 0.05)"
) +
theme(
plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
axis.title = element_text(size = 14),
axis.text = element_text(size = 12),
strip.text = element_text(size = 12),
legend.position = "bottom"
)
return(p)
}
# Example usage
create_volcano_plot(TRNSCRPT_BAT_DA)We provide the following functions to replicate our timewise
and training differential analysis results for each dataset.
Look at the corresponding documentation for each function for more
details, e.g. ?transcript_timewise_da.
proteomics_timewise_da()proteomics_training_da()atac_training_da()atac_timewise_da()immuno_timewise_da()immuno_training_da()metab_timewise_da()metab_training_da()metab_meta_regression()rrbs_differential_analysis()transcript_timewise_da()transcript_training_da()Here, we replicate the protein acetylation timewise differential analysis results provided in MotrpacRatTraining6moData.
timewiseList = list()
for (tissue in c("HEART","LIVER")){
timewiseList[[tissue]] = proteomics_timewise_da("ACETYL", tissue)
}
#> ACETYL_HEART_NORM_DATA
#> Warning: Partial NA coefficients for 157 probe(s)
#> ACETYL_LIVER_NORM_DATA
#> Warning: Partial NA coefficients for 355 probe(s)timewise = do.call("rbind", timewiseList)
# merge with version of results in MotrpacRatTraining6moData
original = combine_da_results(assays="ACETYL")
#> ACETYL_HEART_DA
#> ACETYL_LIVER_DAmerged = merge(original, timewise,
by=c("feature_ID","assay","assay_code","tissue","tissue_code","sex","comparison_group"),
suffixes=c("_orginal","_reproduced"))
# plot
plot(-log10(merged$p_value_orginal),
-log10(merged$p_value_reproduced),
type="p", cex=0.5,
xlab="Original timewise p-value (-log10)",
ylab="Reproduced timewise p-value (-log10)",
main="Timewise DA results for ACETYL data")We can plot the differential analysis results for a single feature
using plot_feature_logfc(). All of the following examples
are different ways to plot results for the same feature.
plot_feature_logfc(assay = "PROT",
tissue = "SKM-GN",
feature_ID = "NP_786937.1",
scale_x_by_time = FALSE)
plot_feature_logfc(assay = "PROT",
tissue = "SKM-GN",
feature_ID = "NP_786937.1",
facet_by_sex = TRUE,
add_adj_p = FALSE)For a detail tutorial on the graphical clustering analysis, please, visit the original vignette part of the MotrpacRatTraining6mo package here.
For questions, bug reporting, and data requests for the
MotrpacRatTraining6moData package, please submit
a new issue here and include as many details as possible.
If the issue is related to functions provided in the
MotrpacRatTraining6mo package, please submit an issue here.
MoTrPAC is supported by the National Institutes of Health (NIH) Common Fund through cooperative agreements managed by the National Institute of Diabetes and Digestive and Kidney Diseases (NIDDK), National Institute of Arthritis and Musculoskeletal Diseases (NIAMS), and National Institute on Aging (NIA).
sessionInfo()
#> R version 4.3.3 Patched (2024-04-09 r86521)
#> Platform: x86_64-apple-darwin20 (64-bit)
#> Running under: macOS Sonoma 14.7.1
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
#>
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> time zone: America/New_York
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] viridis_0.6.5 viridisLite_0.4.2
#> [3] pheatmap_1.0.12 lubridate_1.9.3
#> [5] forcats_1.0.0 stringr_1.5.1
#> [7] purrr_1.0.2 readr_2.1.5
#> [9] tidyr_1.3.1 tibble_3.2.1
#> [11] tidyverse_2.0.0 impute_1.76.0
#> [13] dplyr_1.1.4 ggplot2_3.5.1
#> [15] MotrpacRatTraining6mo_1.6.5 MotrpacRatTraining6moData_2.0.0
#>
#> loaded via a namespace (and not attached):
#> [1] Rdpack_2.6 mnormt_2.1.1 gridExtra_2.3
#> [4] sandwich_3.1-0 rlang_1.1.4 magrittr_2.0.3
#> [7] multcomp_1.4-25 qqconf_1.3.2 compiler_4.3.3
#> [10] vctrs_0.6.5 pkgconfig_2.0.3 crayon_1.5.3
#> [13] fastmap_1.2.0 labeling_0.4.3 ggraph_2.2.1
#> [16] utf8_1.2.4 r2r_0.1.1 rmarkdown_2.27
#> [19] tzdb_0.4.0 xfun_0.45 cachem_1.1.0
#> [22] jsonlite_1.8.8 highr_0.11 tweenr_2.0.3
#> [25] R6_2.5.1 bslib_0.7.0 stringi_1.8.4
#> [28] RColorBrewer_1.1-3 limma_3.58.1 mutoss_0.1-13
#> [31] jquerylib_0.1.4 numDeriv_2016.8-1.1 Rcpp_1.0.12
#> [34] knitr_1.47 zoo_1.8-12 Matrix_1.6-5
#> [37] splines_4.3.3 igraph_2.0.3 timechange_0.3.0
#> [40] tidyselect_1.2.1 rstudioapi_0.16.0 yaml_2.3.8
#> [43] codetools_0.2-20 lattice_0.22-6 Biobase_2.62.0
#> [46] withr_3.0.0 evaluate_0.24.0 survival_3.7-0
#> [49] polyclip_1.10-6 pillar_1.9.0 stats4_4.3.3
#> [52] sn_2.1.1 generics_0.1.3 mathjaxr_1.6-0
#> [55] hms_1.1.3 munsell_0.5.1 scales_1.3.0
#> [58] TFisher_0.2.0 glue_1.7.0 tools_4.3.3
#> [61] metap_1.10 data.table_1.15.4 locfit_1.5-9.10
#> [64] visNetwork_2.1.2 mvtnorm_1.2-5 graphlayouts_1.1.1
#> [67] tidygraph_1.3.1 grid_4.3.3 plotrix_3.8-6
#> [70] rbibutils_2.2.16 edgeR_4.0.16 colorspace_2.1-0
#> [73] ggforce_0.4.2 cli_3.6.3 fansi_1.0.6
#> [76] gtable_0.3.5 sass_0.4.9 digest_0.6.36
#> [79] BiocGenerics_0.48.1 ggrepel_0.9.5 TH.data_1.1-2
#> [82] htmlwidgets_1.6.4 farver_2.1.2 memoise_2.0.1
#> [85] htmltools_0.5.8.1 multtest_2.58.0 lifecycle_1.0.4
#> [88] statmod_1.5.0 MASS_7.3-60.0.1